Renormalization algorithms for Quantum-Many Body Systems in two and higher 
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We describe quantum many-body systems in terms of projected entangled-pair states, which 
naturally extend matrix product states to two and more dimensions. We present an algorithm to 
determine correlation functions in an efficient way. We use this result to build powerful numerical 
simulation techniques to describe the ground state, finite temperature, and evolution of spin systems 
in two and higher dimensions. 

PACS numbers: PACS 



c 

o 
o 



> 

O 

o 

o 

-i— > 



T3 

c 

o 
o 



X 



The theoretical investigation of strongly correlated sys- 
tems is one the most challenging tasks in several fields of 
Physics. Even though several analytical techniques and 
numerical methods have been developed during the last 
decades, there still exist a rich variety of systems which 
remain untractable. Even some of the simplest systems, 
which deal with spins in lattices with short range interac- 
tions, are very hard to simulate numerically. The devel- 
opment of powerful numerical techniques would allow us 
to discover a rich variety of intriguing phenomena and to 
confirm some of the predictions which have been made. 

In the case of ID systems, much analytical insight has 
been gained by finding exact expressions for the ground 
and/or excited eigenstates of some particular Hamiltoni- 
ans, as it is the case for the 1D-AKLT states pj. On 
the other hand, a very powerful numerical simulation 
method known as DMRG Q has allowed us to deter- 
mine physical properties of generic spin chains to an un- 
precedented accuracy. Recent work has also shown how 
DMRG can be adapted to simulate the spin dynamics at 
zero-temperature |3| or at finite temperature and in the 
presence of dissipation 0, [fj . The success of the DMRG 
method and its extensions relies on the existence of the 
so-called matrix product states (MPS) 0, which gener- 
alize the 1D-AKLT states. The DMRG method can be 
understood as a variational method within these MPS 
0, HI 01 , and part of its success relays on the fact that 
correlation functions can be efficiently calculated. 

In two or higher dimensions, however, almost no mod- 
els have been solved exactly. A generalization of DMRG 
to higher dimensions is hard to scale, as the MPS-ansatz 
explicitly assumes a ID configuration. The Monte Carlo 
method 0, on the other hand, is very useful to de- 
scribe certain systems, but for models with frustration 
is, unlike DMRG, plagued by the so-called sign problem. 
The physics of 2D spin systems is therefore not very well 
understood as compared to ID systems; this is very un- 
fortunate as a good understanding would shed new light 
on many open questions in condensed matter, such as 
high-T c superconductivity. 

In this paper, we present a natural generalization of the 
ID MPS to two and higher dimensions and build simu- 



lation techniques based on those states which effectively 
extend DMRG to higher dimensions. We call those states 
projected entangled-pair states (PEPS), since they can 
be understood in terms of pairs of maximally entangled 
states of some auxiliary systems, and that are projected 
in some low-dimensional subspaces locally. This class 
of states includes the generalizations of the 2D AKLT- 
states known as tensor product stat es 11 1 ) which have 
been used for 2D problems (see also |l2l Il3l Kty) but is 
much broader since every state can be represented as a 
PEPS (as long as the dimension of the entangled pairs is 
large enough). We also develop an efficient algorithm to 
calculate correlation functions of these PEPS, and which 
allows us to extend the ID algorithms QEJOIQ to higher 
dimensions. This leads to many interesting applications, 
such as scalable variational methods for finding ground 
or thermal states of spin systems in higher dimensions as 
well as to simulate their time-evolution. For the sake of 
simplicity, we will present our results for a square lattice 
in 2D, but they are easily generalized to higher dimen- 
sions and other geometries. 

Let us start by recalling the representation introduced 
in H of the state VP of TV d-dimensional systems in terms 
of MPS. For that, we substitute each physical system k by 
two auxiliary systems and bk of dimension D (except 
at the extremes of the chain). Systems bk and au+i arc 
in a maximally entangled state \(f>) — X) n =i l n ' n )' wm ch 
is represented in Fig. Ufa) by a solid line (bond) join- 
ing them. The state 'J is obtained by applying a linear 
operator Qk to each pair ak,bk that maps the auxiliary 
systems onto the physical systems, i.e. 
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where the matrices A s k have elements (A s k )i, r = 
(s\Qk\l,r). Note that the indices I and r of each matrix 
A| are related to the left and right bonds of the auxil- 
iary systems with their neighbors, whereas the index s 
denotes the state of the physical system. The function 
Fi is just the trace of the product of the matrices, i.e. it 
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FIG. 1: Graphical representation of MPS in 1 dimension (a), 
in 2 dimensions (b), and of PEPS (c). The bonds represent 
pairs of maximally entangled D-dimensional auxiliary spins 
and the circles denote projectors that map the inner auxiliary 
spins to the physical ones. 



contracts the indices l,r of the matrices A according to 
the bonds. 

As shown in 0, every state can be represented in the 
form JIJ as long as the dimension D can be chosen suf- 
ficiently large. Note, however, that the above picture of 
the state is basically one dimensional, since each auxil- 
iary system is entangled only to one nearest neighbor. 
Thus, these states appear to be better suited to describe 
ID systems, with short range interactions, since a small 
local dimension D may give a good approximation to the 
real state of the whole system. Note also that, as it is 
clear from the above representation, any block of systems 
is only entangled to the rest by at most two maximally 
entangled state of the auxiliary particles and thus its en- 
tropy is bounded by 21og 2 -D, independent of the block 
size. This has been identified as the main reason why 
DMRG does not describe well critical systems, where the 
entropy grows with the logarithm of the block size . 

States in the form Q have also been used to represent 
2D systems 0|. For simplicity let us thus consider a 
2D square lattice of N := Nh x N v systems. The idea 
there is to numerate them in such a way that they can 
be regarded as a long ID system [Fig. db)]. In general, 
this method cannot be extended to larger systems since 
nearest neighbor interactions in the original 2D system 
(for example between 11 and 20) give rise to long inter- 
actions in the effective ID description. Moreover, the 
entropy of some blocks does not scale as the area of the 
block, as it is expected for 2D configurations. For exam- 
ple, the block formed by systems from 6-15 has at most 
a constant entropy of 2 log 2 D. 

For 2D systems we propose to use the description based 
on Fig. Hfc). Each physical system of coordinates (h,v) 
is represented by four auxiliary systems dh, v , bh, v , Ck <v , 
and dh, v of dimension D (except at the borders of the 
lattices). Each of those systems is in a maximally en- 



tangled state 4> with one of its neighbor, as shown in the 
figure. The state ^ is obtained by applying to each site 
one operator Qh,v that maps the auxiliary systems onto 
the physical systems: 
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Here, the A's are four index tensors with elements 
{Af l v ) Uyd j,r = (s\Q h ,v\u,d,l,r). As in the ID case, we 
associate each index of such tensors to each direction (up, 
down, left, and right). Thus, the position with coordi- 
nates (h,v) is represented by a tensor (A s h v ) u ,d,i,r whose 
index s is represents the physical system, and the other 
four indices are associated with the bonds between the 
auxiliary systems and the neighboring ones. The func- 
tion F 2 contracts all these indices u,d,l,r of all tensors 
A according to those bonds. Note that we can general- 
ize this construction to any lattice shape and dimension, 
and that using the construction of |(| one can show that 
any state can be written as a PEPS. In this way, we also 
resolve the problem of the entropy of blocks mentioned 
above, since now this entropy is proportional to the bonds 
that connect such block with the rest, and therefore to 
the area of the block. Note also that, in analogy to the 
MPS 0, the PEPS are guaranteed to be ground states 
of local Hamiltonians. 

We show now how to determine expectation values of 
operators in the state ^ @- We consider a general op- 
erator O — Y[h c Oh.c and define the four-indices tensor 
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where the indices are combined in pairs, i.e., u = 
(u,u'),d = (d,d'),l = (1,1'), and f = (r,r'). One can 
easily show that (*|0|*) = F 2 (E 0h J. Thus, the evalu- 
ation of expectation values consists of contracting indices 
of the tensors E. In order to show how to do this in prac- 
tice, we notice that the tensors associated to the first and 
last rows, once contracted, can be reexpressed in terms 
of a MPS. In particular, we define [compare QJ] 

D 2 

\Ui) := Yl Fi(E d 1 \---E d 1 N N)\di---d N ), (4a) 

D 2 

(U Nv \ := J2 Fi(K\i---E u N : N )(u 1 ...u N \m 

U\ . . .Ujv — 1 

Here we have used the short-hand notation Eh tC := 
Eo h c , and the fact that the tensors in the first and last 
rows have at most three indices [see Fig. HJ:]. Thus, the 
horizontal indices (I, r) of the tensors play the role of the 
indices of each matrix, whereas the vertical ones (d) plays 
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the role of the indices corresponding to the physical sys- 
tems in ID. Analogously, the rows 2, 3, . . . , N v — 1 can be 
considered as matrix product operators (MPO) 

D 2 

U k := J2 F(E^T...,Ef^)\di...d N )( Ul ...u N \. 

di ,ui . . . — 1 

We have (*\0\V) = {U N \U N -i . . . l/ 2 |E7i). 

The evaluation of expectation values poses a serious 
problem since the number of indices proliferate after each 
contraction. For example, the vector \U<2) := U2\U\) can 
be written as the MPS 14a|) but with the substitution 

D 2 

E t^ T. E t® E i^ n - ( 5 ) 

d»=l 

This last tensor has more (right and left) indices than the 
original one. Thus, every time we apply the MPO Uk to 
a MPS \Uk-i) the number of indices increases, and thus 
the problem soon becomes intractable. Now we introduce 
a numerical algorithm inspired by DMRG to numerically 
determine F2 {Eo h c ) and to overcome this problem. 

Given an unnormalized MPS \^>a) parameterized by 
D x D matrices {^j, the goal is to find another MPS 
\i/>b)> parameterized by Df x Df matrices {-Bf}, where 
Df < D is a prescribed number. This has to be done 
such that K := \\\ipA) ~ l^s)l| 2 1S minimal, i.e., such 
that that |V>s) gives the best approximation to \4>a)- We 
have developed an algorithm that achieves this task in 
an iterative way. The key insight is that K is quadratic 
in all components of the matrices {Bj^}, and hence if all 
these matrices are fixed except one of them (say 73* ) K 
is quadratic in the components of B? ; the optimal choice 
for B? thus amounts to solving a trivial system of linear 
equations. Having done that, one moves to the next site 
i + 1, fixes all other ones and repeats the same procedure. 
After a few sweeps back and forth the optimal MPS is 
found. Note that the error in the approximation is ex- 
actly known and if it becomes too large one can always 
increase Df ; in all relevant situations we encountered the 
error could be made very small even with moderate Df. 
The same reasoning holds for MPS defined with periodic 
instead of open boundary conditions. In this latter case 
considered here, one can further simplify the system of 
linear equations by performing a singular value decompo- 
sition of Bj and keeping only one of the unitary matrices 
at each step, analogously as one does in DMRG. 

Thus, in order to evaluate an arbitrary expectation 
value we first determine the MPS If/2) which is the clos- 
est to U2\Ui) but with a fixed dimensions Df of the cor- 
responding matrices. Then, we determine If/3), which 
is the closest to C/3 IC/2) 5 and continue in this vein until 
we finally determine (*|0|*) ~ (U N \U N -i). Interest- 
ingly enough, this method to calculate expectation val- 
ues and to determine optimal approximations to MPS 



can be adapted to develop very efficient algorithms to 
determine the ground states of 2D Hamiltonians and the 
time evolution of PEPS by extending DMRG and the 
time evolution schemes to 2D. 

Let us start with an algorithm to determine the ground 
state of a Hamiltonian with short range interactions on 
a square Nh x N v lattice. The goal is to determine the 
PEPS (0 with a given dimension D which minimizes the 
energy. Following the idea is to iteratively optimize 
the tensors A a h c one by one while fixing all the other ones. 
The crucial observation is the fact that the exact energy 
of \if>) (and also its normalization) is a quadratic function 
of the components of the tensor A| _ to be optimized, 
which we write as a vector x; hence the energy can be 
expressed in terms of an effective Hamiltonian: 

E = ^1" a (6) 
x^Nx 

The denominator takes the normalization of the state 
into account. This expression can readily be minimized 
as it is equivalent to a generalized eigenvalue problem. 
It turns out that H e ff and N can be efficiently evalu- 
ated by the methods described above. In the case of N, 
the MPS |f/i) (ptaj) constructed from E-^ i can be prop- 
agated up to the row v — 1 with the technique outlined 
before. Similarly, the last row (f/jvj <|4b(l can be prop- 
agated up to row v + 1. The tensors E-^ can now be 
contracted with these two MPS from h = l..h — 1, and 
similarly from h — N^.-h + 1. The remaining tensor has 
4 (double) indices from which one can readily determine 
N. H e ff can be determined in an analogous way, but 
here the procedure has to be repeated for every term in 
the Hamiltonian (i.e. in the order of 2Ny l N v times in the 
case of nearest neighbor interactions) . Thus both TV and 
H e f f can be calculated efficiently. Therefore the optimal 
A%- can be determined, and one can proceed with the 
following projector, iterating the procedure until conver- 
gence. 

In a practical implementation one can save much time 
by storing appropriate tensors, implementing the algo- 
rithm in a parallel way, doing sparse tensor multiplica- 
tions, and making use of quantum numbers and reflec- 
tion symmetries. A more efficient variant can also be 
constructed by an iterative procedure which resembles 
the infinite-size DMRG-algorithm, where new rows are 
inserted in the middle of the lattice. 

Let us next move to describe how time-evolution can be 
simulated on a PEPS. We will assume that the Hamilto- 
nian only couples nearest neighbors, although more gen- 
eral settings can be considered. The simplest scheme 
would work by optimally mapping a given PEPS to an- 
other PEPS after an infinitesimal time-step 1 — iHSt. It 
can readily be checked that, up to first order of 6t, the 
action of this operator is to map a TPS with auxiliary 
dimension D onto a new TPS with dimension of the aux- 
iliary bonds nD\ here n represents the minimal number 
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of terms needed to express the couplings as tensor prod- 
ucts of local operators plus 1 (e.g. n = 2 for the Ising 
interaction and n = 4 for the Heisenberg interaction) . In 
analogy to the method introduced above for MPS, one 
can approximate this new PEPS with another one hav- 
ing again bonds of auxiliary dimension D. The algorithm 
to achieve this is a direct generalization of the method 
introduced to reduce the D of MPS: again several sweeps 
over all projectors have to be done, and the only differ- 
ence is that at each step correlation functions of a PEPS 
have to be calculated instead of correlations function of 
a MPS. This can be done using the methods introduced 
before. Of course there are again many possibilities to 
boost the accuracy and to reduce the computational cost 
of such an implementation, such as using the Trotter de- 
composition as in Q and then using the sweeps to opti- 
mize the state. This algorithm can also be used to solve 
finite temperature or dissipation problems by extending 
the ideas of and 0. 

Let us now illustrate our methods with an example. 
We consider a 2D lattice of spin 1/2-particles where near- 
est neighbors interact via the antiferromagnetic Heisen- 
berg interaction with coupling constant J = 1. We use 
the time-evolution algorithm for evolving the PEPS in 
imaginary time; in this way we illustrate both the fact 
that the new formalism allows us to find ground states as 
well as to describe time-evolution. We implemented the 
algorithm as follows: we start with a completely separa- 
ble state |^o) m which the spins are rotated by an angle 
7r/16 with respect to the previous one, and which can 
trivially be written as a PEPS. Using the Trotter decom- 
position, we divide each time step into 4 parts in which 
each spin is only interacting with one neighbor; as we 
are considering the Heisenberg interaction, the dimen- 
sion D between the 2 interacting spins gets multiplied 
by a factor of 4. Let us parameterize this new PEPS 
with the corresponding tensors Bf t v . After each of these 
substeps, we want to reduce the dimension again to the 
original one giving rise to the PEPS that optimally 
approximates the exact B s h . This is done in an iterative 
way, row by row, until convergence. Fixing all rows but 
one, the problem of finding the optimal projectors in this 
row is equivalent to the problem of approximating a MPS 
with another one with lower dimension (the physical di- 
mension of the MPS is the product of the bonds going up 
and down) , which can on itself done in an iterative way 
as outlined above. Note that the computational cost of 
the algorithm is polynomial in N and D. 

We have first considered a 4 x 4 lattice on which the 
imaginary time evolution can be determined exactly. In 
Fig. 2, we plotted the exact evolution versus the one 
where the evolution is approximated variationally within 
the PEPS with bonds of dimension D = 2, 3 (D = 4 can- 
not be distinguished from the exact result). We used the 
same Trotter approximation for the exact and variational 
simulations with St = — 3i/100. It is remarkable that 
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FIG. 2: Imaginary time evolution with the Heisenberg and 
a frustrated Heisenberg interaction on a 4 x 4 lattice, and 
D = 2,Df = 16 (dotted) and D = 3, D f = 35 (dashed); 
the D = 3 results are almost indistinguishable from the exact 
ones (full line). The insert presents the evolution for D — 2 
on a 10 x 10 lattice. 



even for D = 3 we obtain a very good approximations, 
both regarding time evolution and ground state energy. 
The algorithm clearly converges to the ground state, and 
the difference between the exact ground state energy and 
the one obtained with our scalable algorithm rapidly de- 
creases with L>[l3; more specifically, 1 — E var /E exact is 
given by .35, .02; .004; 0.0008 for D = 1, 2; 3; 4 (note that 
the trivial situation D = 1 corresponds to the Neel state) . 
We also repeated the same simulation but with a frus- 
trated Heisenberg Hamiltonian, obtained by making 1 
out of every 4 interactions on each spin ferromagnetic in- 
stead of antiferromagnetic. Again very good agreement 
with the exact results is obtained; note that the energy 
converges more slowly due to the fact that the energy gap 
is smaller. The insert of Fig. @ presents some simula- 
tion results for the imaginary time evolution for a square 
10 x 10 lattice for both the Heisenberg antiferromagnet 
and the frustrated case. The convergence is again very 
fast, and increasing D from 2 to 3 (not shown in the 
plot) allows us to find a better value for the energy of 
the ground state. Note that we can easily handle larger 
systems and, using the appropriate numerical techniques, 
eventually increase the value of D. 

In conclusion, we have introduced the class of PEPS 
and showed how they arise naturally in the context of 
constructing variational ground states for spin Hamilto- 
nians on higher dimensional lattices. We presented an 
efficient algorithm for calculating correlation functions, 
which leads to scalable variational methods for finding 
ground states and for describing their real or imaginary 
time evolution. Interestingly, the methods described also 
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apply in the case of different geometries, of evolution in 
the presence of dissipation, and for finding finite-T states. 
It is also possible to identify quite generic classes of PEPS 
for which 2-point correlation functions can be calculated 
analytically QJ3. We also note that the concept of PEPS 
could be very useful for the description of 2-dimensional 
transport problems, as the PEPS generalize the matrix 
product states which proved to be very useful in the 1-D 
case [l^ . 
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